{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "# utility modules\n",
    "import glob\n",
    "import os\n",
    "import sys\n",
    "import re\n",
    "import io\n",
    "\n",
    "# the usual suspects\n",
    "import numpy as np\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.io import shapereader\n",
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "import cartopy.io.img_tiles as cimgt\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "# import matplotlib.image as mpimg\n",
    "\n",
    "# modules you'll need if you're downloading the data\n",
    "from icepyx import icesat2data as ipd\n",
    "import shutil\n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "from shapely.geometry import Point, Polygon\n",
    "\n",
    "import fiona\n",
    "import pyproj\n",
    "import h5py\n",
    "\n",
    "# run matplotlib in 'widget' mode\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "53.09700798273372\n"
     ]
    }
   ],
   "source": [
    "from readers.read_HDF5_ATL06 import atl06_to_dict\n",
    "\n",
    "ATL06_file = 'ATL06_20190729075244_04770406_003_01.h5'\n",
    "beam = ['gt1l','gt1r','gt2l','gt2r','gt3l','gt3r']\n",
    "\n",
    "# load kml data to select points inside glacier area\n",
    "fiona.drvsupport.supported_drivers['kml'] = 'rw' \n",
    "fiona.drvsupport.supported_drivers['KML'] = 'rw' \n",
    "glacier_filepath = os.getcwd()+'/sikkimglaciers.kml'\n",
    "gdfg = gpd.read_file(glacier_filepath) \n",
    "mask = gdfg.geometry.unary_union\n",
    "\n",
    "D6 = {} \n",
    "thickness = 0\n",
    "for i in beam:\n",
    "    D6[i] = {}\n",
    "    data = atl06_to_dict(ATL06_file, i, index=None, epsg=3031)\n",
    "    flag = np.zeros(data['longitude'].size, dtype=bool)\n",
    "    for j in range(0,data['longitude'].size):\n",
    "        flag[j] = Point(data['longitude'][j],data['latitude'][j]).within(mask) \n",
    "    D6[i]['lon'] = data['longitude'][flag]\n",
    "    D6[i]['lat'] = data['latitude'][flag]\n",
    "    D6[i]['h_ortho'] = data['h_li'][flag]+data['geoid_h'][flag]\n",
    "    D6[i]['thickness'] = data['h_li'][flag]+data['geoid_h'][flag]-data['dem_h'][flag]\n",
    "    # calculate the mean glacier thickness\n",
    "    ext = (D6[i]['thickness'][:]>0) & (D6[i]['thickness'][:]<1000)\n",
    "    D6[i]['filter_thickness'] = D6[i]['thickness'][ext]\n",
    "    thickness = thickness + np.average(D6[i]['filter_thickness'][~np.isnan(D6[i]['filter_thickness'])])\n",
    "    # standard deviation gives no meaning because the glacier thickness varies greatly from the edge to center\n",
    "thickness = thickness/len(beam)\n",
    "print(thickness)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fd884634b92a444b8b121821549070a5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the points\n",
    "extent = [88.1, 88.35, 27.71, 28.1]\n",
    "request = cimgt.OSM()    \n",
    "fig = plt.figure(figsize=(9, 13))\n",
    "ax = plt.axes(projection=request.crs)\n",
    "cmap = cm.viridis_r\n",
    "cmap.set_bad(alpha=0.0)\n",
    "for gtx in ['gt1l','gt1r','gt2l','gt2r','gt3l','gt3r']:\n",
    "    flag = (D6[gtx]['thickness'][:]>0) & (D6[gtx]['thickness'][:]<1000)\n",
    "    sc = ax.scatter(D6[gtx]['lon'][flag],D6[gtx]['lat'][flag],\n",
    "        c=D6[gtx]['thickness'][flag], s=5, cmap=cmap, transform=ccrs.PlateCarree())  \n",
    "cbar = fig.colorbar(sc,ax=ax,extend='both',extendfrac=0.0375, pad=0.03, drawedges=False)\n",
    "cbar.solids.set_rasterized(True)\n",
    "cbar.ax.tick_params(which='both', length=16, width=1, direction='in')\n",
    "cbar.ax.set_ylabel('Land ice thickness')\n",
    "cbar.ax.set_xlabel('m')\n",
    "cbar.ax.xaxis.set_label_coords(0.50,1.04)      \n",
    "gl = ax.gridlines(draw_labels=True, alpha=0.2)\n",
    "# gl.xlabels_top = gl.ylabels_right = False\n",
    "gl.xformatter = LONGITUDE_FORMATTER\n",
    "gl.yformatter = LATITUDE_FORMATTER\n",
    "ax.set_extent(extent)\n",
    "ax.add_image(request, 10)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "for gtx in ['gt1l','gt1r','gt2l','gt2r','gt3l','gt3r']:\n",
    "    flag = (D6[gtx]['thickness'][:]>0) & (D6[gtx]['thickness'][:]<1000)\n",
    "    y = plt.scatter(D6[gtx]['lat'][flag],D6[gtx]['thickness'][flag])\n",
    "plt.show()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
